//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007, 2008 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <assert.h>

#include "ephi.hpp"

#define EMOD (SPEED_OF_LIGHT*SPEED_OF_LIGHT*MAGNETIC_CONSTANT)
#define BMOD MAGNETIC_CONSTANT

PMap::PMap (const vect3d& dim, prec_t resolution, prec_t contact, prec_t smooth, size_t smooth_count)
{
  xc = 2 * prec2size_t(dim.x / resolution);
  yc = 2 * prec2size_t(dim.y / resolution);
  zc = 2 * prec2size_t(dim.z / resolution);
  hash.resize (prec2size_t (xc * yc * zc * 1.33), 0);
  bs = vect3d (resolution, resolution, resolution);
  hbs = 0.5 * bs;
  rbs = bs.recip();
  corner = vect3d (-(xc/2 * bs.x), -(yc/2 * bs.y), -(zc/2 * bs.z));
  cmag = contact * contact;
  acc = (5 * resolution) * (5 * resolution);
  //acc = (3 * bs).magnitude();
  this->smooth = smooth;
  this->smooth_count = smooth_count;
  fsmooth = 0.1 * smooth;
  smag = smooth * smooth;

  xs = prec2size_t(ceil(smooth / resolution)) * 2 + 1;
  ys = prec2size_t(ceil(smooth / resolution)) * 2 + 1;
  zs = prec2size_t(ceil(smooth / resolution)) * 2 + 1;
  xss = xs / 2;
  yss = ys / 2;
  zss = zs / 2;

  prec_t sf = smooth / 8.0;
  for ( size_t i = 0; i < 8; i++ )
    {
      sv[i][0] = (i + 1) * sf;
      sv[i][1] = sv[i][0] * sv[i][0];
    }

  prec_t R = smooth;
  prec_t R2 = R * R;
  prec_t R3 = R2 * R;
  prec_t R5 = R3 * R2;
  prec_t R7 = R5 * R2;

  // M
  f35_8R3 = 35 / (8 * R3);
  f21_4R5 = 21 / (4 * R5);
  f15_8R7 = 15 / (8 * R7);
  // N
  f21_2R5 = 21 / (2 * R5);
  f15_2R7 = 15 / (2 * R7);
  // K
  f15_R7 = 15 / R7;
  // pot
  f35_16R = 35 / (16 * R);
  f35_16R3 = 35 / (16 * R3);
  f21_16R5 = 21 / (16 * R5);
  f5_16R7 = 5 / (16 * R7);
}

PMap::~PMap ()
{
}

bool PMap::addQ (const vect3d& pos, const vect3d& v, prec_t q)
{
  vect3d n = pos - corner;
  vect3d tp = vect3d (n.x * rbs.x, n.y * rbs.y, n.z * rbs.z);
  vect3d tpf = vect3d (floor(tp.x), floor(tp.y), floor(tp.z));
  size_t xi = prec2size_t(tpf.x);
  size_t yi = prec2size_t(tpf.y);
  size_t zi = prec2size_t(tpf.z);
  if (! (xi < xc && yi < yc && zi < zc) )
      return false;
  size_t index = zi + yi * zc + xi * yc * zc;
  size_t i = index % hash.size();
  cell_t *cell = findCell (index);
  if ( !cell )
    {
      cell = a_cells.getz();
      cell->mid = corner + vect3d (xi * bs.x, yi * bs.y, zi * bs.z) + hbs;
      cell->index = index;
      cell->next = hash[i];
      hash[i] = cell;
    }
  pnode_t *pn = a_pnodes.get();
  pn->p.pos = pos;
  pn->p.v = v;
  pn->p.q = q;
  pn->p.dobf = false;
  if ( v.x != 0 || v.y != 0 || v.z != 0 )
    {
      pn->p.dobf = true;
      cell->dobf = true;
    }
  cell->add (pn);
  return true;
}

void PMap::addUniformSphere (const vect3d& pos, prec_t r, prec_t q)
{
  spheres.push_back (sphere_t());
  spheres.back().pos = pos;
  spheres.back().r = r;
  spheres.back().rr = r * r;
  spheres.back().q = q;
  spheres.back().q_by_r = q / r;
  spheres.back().q_by_r3 = q / (r * r * r);
  spheres.back().rho = 3 * q / (4.0 * PREC_PI * r * r * r);
}

bool PMap::inContact (const vect3d& pos)
{
  if ( cmag == 0 )
      return false;
  vect3d n = pos - corner;
  vect3d tp = vect3d (n.x * rbs.x, n.y * rbs.y, n.z * rbs.z);
  vect3d tpf = vect3d (floor(tp.x), floor(tp.y), floor(tp.z));
  size_t xi = prec2size_t(tpf.x);
  size_t yi = prec2size_t(tpf.y);
  size_t zi = prec2size_t(tpf.z);
  for ( size_t i = 0; i < 3; i++ )
    {
      for ( size_t j = 0; j < 3; j++ )
        {
          for ( size_t k = 0; k < 3; k++ )
            {
              cell_t *cell = findCell(makeIndex (xi + i - 1, yi + j - 1, zi + k - 1));
              if ( !cell )
                  continue;
              pnode_t *pn = cell->pnodes;
              while ( pn )
                {
                  if ( (pn->p.pos - pos).magnitude() < cmag )
                      return true;
                  pn = pn->next;
                }
            }
        }
    }
  return false;
}

void PMap::getPos (const vect3d& pos, size_t& xi, size_t& yi, size_t& zi)
{
  vect3d n = pos - corner;
  vect3d tp = vect3d (n.x * rbs.x, n.y * rbs.y, n.z * rbs.z);
  vect3d tpf = vect3d (floor(tp.x), floor(tp.y), floor(tp.z));
  xi = prec2size_t(tpf.x);
  yi = prec2size_t(tpf.y);
  zi = prec2size_t(tpf.z);
}

PMap::cell_t * PMap::findCell (size_t index)
{
  size_t i = index % hash.size();
  cell_t *cell = hash[i];
  while ( cell && cell->index != index )
      cell = cell->next;
  return cell;
}

void PMap::recalc ()
{
  size_t cc = 0;
  for ( size_t i = 0; i < hash.size(); i++ )
    {
      cell_t * cell = hash[i];
      while ( cell )
        {
          cc++;
          cell->qsum = 0;
          cell->qmid.clear();
          cell->pcount = 0;
          pnode_t *pn = cell->pnodes;
          while ( pn )
            {
              cell->qsum += pn->p.q;
              cell->j += pn->p.q * pn->p.v;
              cell->pcount++;
              pn = pn->next;
            }
          cell->particles = a_particles.get(cell->pcount);
          cell->jmid.clear();
          prec_t jadd = 0;
          pn = cell->pnodes;
          size_t pi = 0;
          while ( pn )
            {
              cell->particles[pi] = pn->p;
              if ( cell->qsum != 0 )
                  cell->qmid += (pn->p.q / cell->qsum) * pn->p.pos;
              prec_t jf = pn->p.q * (pn->p.v * cell->j);
              if ( jf > 0 )
                {
                  cell->jmid += jf * pn->p.pos;
                  jadd += jf;
                }
              pn = pn->next;
              pi++;
            }
          if ( jadd != 0 )
              cell->jmid /= jadd;
          else
              cell->jmid = cell->qmid;
          cell = cell->next;
        }
    }

  cells.reserve (cc);

  for ( size_t i = 0; i < hash.size(); i++ )
    {
      cell_t * cell = hash[i];
      while ( cell )
        {
          cells.push_back (cell);
          cell = cell->next;
        }
    }
}

prec_t PMap::getPotential (const vect3d& pos)
{
  CubicScalar epot;
  CubicVector cbf, cef;
  epot.clear();
  addCubicValues (pos, epot, cbf, cef);
  return epot.f;
}

prec_t PMap::getCrudePotential (const vect3d& pos)
{
  prec_t ret = 0;
  for ( size_t i = 0; i < cells.size(); i++ )
      ret += cells[i]->qsum / (pos - cells[i]->qmid).length();
  ret *= EMOD;
  return ret;
}

void PMap::addCrudeEField (const vect3d& pos, vect3d& ef)
{
  vect3d add;
  add.clear();
  for ( size_t i = 0; i < cells.size(); i++ )
    {
      vect3d r =  (pos - cells[i]->qmid);
      prec_t r2 = r * r;
      prec_t rlen = prec_t_sqrt (r2);
      add += (cells[i]->qsum / (r2 * rlen)) * r;
    }
  ef += EMOD * add;
}

prec_t PMap::getChargeDensity (const vect3d& pos)
{
  prec_t ret = 0;
  prec_t mm = (smooth + bs.length());
  mm = mm * mm;
  size_t xi, yi, zi;
  getPos (pos, xi, yi, zi);
  size_t n = 2 * (size_t)ceil(smooth / bs.length()) + 3;
  size_t hn = n / 2;
  for ( size_t i = 0; i < n; i++ )
    {
      for ( size_t j = 0; j < n; j++ )
        {
          for ( size_t k = 0; k < n; k++ )
            {
              cell_t *cell = findCell(makeIndex (xi + i - hn, yi + j - hn, zi + k - hn));
              if ( !cell )
                  continue;
              if ( (cell->mid - pos).magnitude() >= mm )
                  continue;
              for ( size_t l = 0; l < cell->pcount; l++ )
                {
                  particle_t& p = cell->particles[l];
                  prec_t r2 = (pos - p.pos).magnitude();
                  if ( r2 >= smag )
                      continue;
                  prec_t k = (p.q * 105) / (4 * PREC_PI * 8 * smooth * smag);
                  prec_t add = 1.0 - r2 / smag;
                  ret += k * add * add;
                }
            }
        }
    }
  return ret;
}

void PMap::addCubicValues (const vect3d& pos, CubicScalar& out_epot, CubicVector& out_cbf, CubicVector& out_cef)
{
  CubicScalar epot;
  CubicVector cef, cbf;
  epot.clear();
  cef.clear();
  cbf.clear();

#if 0
  // at first calculate the histogram of how many particles are nearby,
  // based on that we can select the smooth factor
  size_t xi, yi, zi;
  getPos (pos, xi, yi, zi);
  size_t hist[7] = {0, 0, 0, 0, 0, 0, 0};
  for ( size_t i = 0; i < xs; i++ )
    {
      for ( size_t j = 0; j < ys; j++ )
        {
          for ( size_t k = 0; k < zs; k++ )
            {
              cell_t *cell = findCell(makeIndex (xi + i - xss, yi + j - yss, zi + k - zss));
              if ( !cell )
                  continue;
              for ( size_t l = 0; l < cell->pcount; l++ )
                {
                  particle_t& p = cell->particles[l];
                  prec_t r2 = (pos - p.pos).magnitude();
                  for ( size_t m = 0; m < 7; m++ )
                    {
                      if ( r2 < sv[m][1] )
                        {
                          hist[m]++;
                          break;
                        }
                    }
                }
            }
        }
    }

  size_t c = 0;
  size_t si;
  for ( si = 0; si < 7; si++ )
    {
      c += hist[si];
      //if ( c != 0 ) printf ("%zu: %zu %zu\n", si, c, hist[si]);
      if ( c > smooth_count )
          break;
    }
  //if ( si != 7 ) printf ("si=%zu, c=%zu\n", si, c);
  si = 7;
  prec_t smooth = sv[si][0];
  prec_t smag = sv[si][1];
#endif

  for ( size_t i = 0; i < cells.size(); i++ )
    {
      cell_t * cell = cells[i];
      vect3d r = pos - cell->qmid;
      prec_t r2 = r.magnitude();
      vect3d jr = pos - cell->jmid;
      prec_t jr2 = jr.magnitude();
      if ( r2 > acc && jr2 > acc )
        {
          prec_t rlen = prec_t_sqrt (r2);
          prec_t r3 = r2 * rlen;
          prec_t r5 = r3 * r2;
          prec_t r7 = r5 * r2;
          prec_t r9 = r7 * r2;

          epot.f += cell->qsum / rlen;
          vect3d ef = (cell->qsum / r3) * r;
          cef.f += ef;
          epot.dfdx += -ef.x;
          epot.dfdy += -ef.y;
          epot.dfdz += -ef.z;
          prec_t f = cell->qsum / r5;
          vect3d dedx = f * (vect3d(r2, 0, 0) - 3 * r.x * r);
          vect3d dedy = f * (vect3d(0, r2, 0) - 3 * r.y * r);
          vect3d dedz = f * (vect3d(0, 0, r2) - 3 * r.z * r);
          cef.dfdx += dedx;
          cef.dfdy += dedy;
          cef.dfdz += dedz;
          epot.d2fdxdy += -dedy.x;
          epot.d2fdxdz += -dedy.x;
          epot.d2fdydz += -dedz.y;
          f = 3 * cell->qsum / r7;
          vect3d d2edxdy = f * (5 * r.x * r.y * r - vect3d (r2 * r.y, r2 * r.x, 0));
          vect3d d2edxdz = f * (5 * r.x * r.z * r - vect3d (r2 * r.z, 0, r2 * r.x));
          vect3d d2edydz = f * (5 * r.y * r.z * r - vect3d (0, r2 * r.z, r2 * r.y));
          cef.d2fdxdy += d2edxdy;
          cef.d2fdxdz += d2edxdz;
          cef.d2fdydz += d2edydz;
          epot.d3fdxdydz += -d2edydz.x;
          f = 15 * cell->qsum / r9;
          vect3d d3edxdydz = f * (vect3d(r.y * r.z * r2, r.x * r.z * r2, r.x * r.y * r2) - 7 * r.x * r.y * r.z * r);
          cef.d3fdxdydz += d3edxdydz;

          if ( cell->dobf )
          {
            prec_t jrlen = prec_t_sqrt (jr2);
            prec_t jr3 = jr2 * jrlen;
            prec_t jr5 = jr3 * jr2;
            prec_t jr7 = jr5 * jr2;
            prec_t jr9 = jr7 * jr2;
            vect3d bf = (1 / jr3) * jr;
            cbf.f += cell->j % bf;

            prec_t f = 1.0 / jr5;
            vect3d dbdx = f * (vect3d(jr2, 0, 0) - 3 * jr.x * jr);
            vect3d dbdy = f * (vect3d(0, jr2, 0) - 3 * jr.y * jr);
            vect3d dbdz = f * (vect3d(0, 0, jr2) - 3 * jr.z * jr);
            cbf.dfdx += cell->j % dbdx;
            cbf.dfdy += cell->j % dbdy;
            cbf.dfdz += cell->j % dbdz;

            f = 3 / jr7;
            vect3d d2bdxdy = f * (5 * jr.x * jr.y * r - vect3d (jr2 * jr.y, jr2 * jr.x, 0));
            vect3d d2bdxdz = f * (5 * jr.x * jr.z * r - vect3d (jr2 * jr.z, 0, jr2 * jr.x));
            vect3d d2bdydz = f * (5 * jr.y * jr.z * r - vect3d (0, jr2 * jr.z, jr2 * jr.y));
            cbf.d2fdxdy += cell->j % d2bdxdy;
            cbf.d2fdxdz += cell->j % d2bdxdz;
            cbf.d2fdydz += cell->j % d2bdydz;
            f = 15 / jr9;
            vect3d d3bdxdydz = f * (vect3d(jr.y * jr.z * jr2, jr.x * jr.z * jr2, jr.x * jr.y * jr2)
                - 7 * jr.x * jr.y * jr.z * r);
            cbf.d3fdxdydz += cell->j % d3bdxdydz;
          }
        }
      else
        {
          for ( size_t j = 0; j < cell->pcount; j++ )
            {
              particle_t& p = cell->particles[j];
              vect3d r = pos - p.pos;
              prec_t r2 = r * r;
              prec_t rlen = prec_t_sqrt (r2);
              bool do_derivs = true;
              if ( rlen < 1e-13 )
                  do_derivs = false;
              bool do_ef = true, do_bf = p.dobf;

              if ( r2 < smag )
                {
                  prec_t r4 = r2 * r2;
                  prec_t r6 = r4 * r2;
                  prec_t M = f35_8R3 - f21_4R5 * r2 + f15_8R7 * r4;
                  vect3d ef = p.q * M * r;
                  prec_t N = -f21_2R5 + f15_2R7 * r2;
                  vect3d dedx = p.q * (r.x * N * r + vect3d(M, 0, 0));
                  vect3d dedy = p.q * (r.y * N * r + vect3d(0, M, 0));
                  vect3d dedz = p.q * (r.z * N * r + vect3d(0, 0, M));
                  prec_t K = -f15_R7;
                  vect3d d2edxdy = p.q * (r.x * r.y * K * r + vect3d (r.y * N, r.x * N, 0));
                  vect3d d2edxdz = p.q * (r.x * r.z * K * r + vect3d (r.z * N, 0, r.x * N));
                  vect3d d2edydz = p.q * (r.y * r.z * K * r + vect3d (0, r.z * N, r.y * N));
                  vect3d d3edxdydz = p.q * (vect3d(r.y * r.z * K, r.x * r.z * K, r.x * r.y * K));

                  cef.f += ef;
                  if ( do_derivs )
                    {
                      cef.dfdx += dedx;
                      cef.dfdy += dedy;
                      cef.dfdz += dedz;
                      cef.d2fdxdy += d2edxdy;
                      cef.d2fdxdz += d2edxdz;
                      cef.d2fdydz += d2edydz;
                      cef.d3fdxdydz += d3edxdydz;
                    }

                  epot.f += p.q * (f35_16R - f35_16R3*r2 + f21_16R5*r4 - f5_16R7*r6);
                  if ( do_derivs )
                    {
                      epot.dfdx += -p.q * r.x * M;
                      epot.dfdy += -p.q * r.y * M;
                      epot.dfdz += -p.q * r.z * M;
                    }

                  do_ef = false;

                  if ( do_bf )
                    {
                      cbf.f += p.v % ef;
                      if ( do_derivs )
                        {
                          cbf.dfdx += p.v % dedx;
                          cbf.dfdy += p.v % dedy;
                          cbf.dfdz += p.v % dedz;
                          cbf.d2fdxdy += p.v % d2edxdy;
                          cbf.d2fdxdz += p.v % d2edxdz;
                          cbf.d2fdydz += p.v % d2edydz;
                          cbf.d3fdxdydz += p.v % d3edxdydz;
                        }
                      do_bf = false;
                    }
                  
#if 0
                  prec_t C = p.q / smag;
                  vect3d ef = (C / rlen) * r;
                  prec_t f = C / r3;
                  vect3d dedx = f * (vect3d(r2, 0 ,0) - r.x * r);
                  vect3d dedy = f * (vect3d(0, r2, 0) - r.y * r);
                  vect3d dedz = f * (vect3d(0, 0, r2) - r.z * r);
                  f = C / r5;
                  vect3d d2edxdy = f * (3 * r.x * r.y * r - vect3d(r.y * r2, r.x * r2, 0));
                  vect3d d2edxdz = f * (3 * r.x * r.z * r - vect3d(r.z * r2, 0, r.x * r2));
                  vect3d d2edydz = f * (3 * r.y * r.z * r - vect3d(0, r.z * r2, r.y * r2));
                  f = 3 * C / r7;
                  vect3d d3edxdydz = f * (vect3d(r.y * r.z * r2, r.x * r.z * r2, r.x * r.y * r2) - 5 * r.x * r.y * r.z * r);

                  cef.f += ef;
                  if ( rlen > fsmooth ) {
                      cef.dfdx += dedx;
                      cef.dfdy += dedy;
                      cef.dfdz += dedz;
                      cef.d2fdxdy += d2edxdy;
                      cef.d2fdxdz += d2edxdz;
                      cef.d2fdydz += d2edydz;
                      cef.d3fdxdydz += d3edxdydz;
                    }

                  epot.f += p.q / smooth; // note: all pot derivatives are 0 here

                  do_ef = false;

                  if ( do_bf )
                    {
                      cbf.f += p.v % ef;
                      if ( rlen > fsmooth ) {
                          cbf.dfdx += p.v % dedx;
                          cbf.dfdy += p.v % dedy;
                          cbf.dfdz += p.v % dedz;
                          cbf.d2fdxdy += p.v % d2edxdy;
                          cbf.d2fdxdz += p.v % d2edxdz;
                          cbf.d2fdydz += p.v % d2edydz;
                          cbf.d3fdxdydz += p.v % d3edxdydz;
                        }
                      do_bf = false;
                    }
#endif
                }

              if ( do_ef || do_bf )
                {
                  prec_t r3 = r2 * rlen;
                  prec_t r5 = r3 * r2;
                  prec_t r7 = r5 * r2;
                  prec_t r9 = r7 * r2;

                  vect3d ef = (p.q / r3) * r;
                  prec_t f = p.q / r5;
                  vect3d dedx = f * (vect3d(r2, 0, 0) - 3 * r.x * r);
                  vect3d dedy = f * (vect3d(0, r2, 0) - 3 * r.y * r);
                  vect3d dedz = f * (vect3d(0, 0, r2) - 3 * r.z * r);
                  f = 3 * p.q / r7;
                  vect3d d2edxdy = f * (5 * r.x * r.y * r - vect3d (r2 * r.y, r2 * r.x, 0));
                  vect3d d2edxdz = f * (5 * r.x * r.z * r - vect3d (r2 * r.z, 0, r2 * r.x));
                  vect3d d2edydz = f * (5 * r.y * r.z * r - vect3d (0, r2 * r.z, r2 * r.y));
                  f = 15 * p.q / r9;
                  vect3d d3edxdydz = f * (vect3d(r.y * r.z * r2, r.x * r.z * r2, r.x * r.y * r2) - 7 * r.x * r.y * r.z * r);

                  if ( do_ef )
                    {
                      cef.f += ef;
                      cef.dfdx += dedx;
                      cef.dfdy += dedy;
                      cef.dfdz += dedz;
                      cef.d2fdxdy += d2edxdy;
                      cef.d2fdxdz += d2edxdz;
                      cef.d2fdydz += d2edydz;
                      cef.d3fdxdydz += d3edxdydz;

                      epot.f += p.q / rlen;
                      epot.dfdx += -ef.x;
                      epot.dfdy += -ef.y;
                      epot.dfdz += -ef.z;
                      epot.d2fdxdy += -dedy.x;
                      epot.d2fdxdz += -dedy.x;
                      epot.d2fdydz += -dedz.y;
                      epot.d3fdxdydz += -d2edydz.x;
                    }

                  if ( do_bf )
                    {
                      cbf.f += p.v % ef;
                      cbf.dfdx += p.v % dedx;
                      cbf.dfdy += p.v % dedy;
                      cbf.dfdz += p.v % dedz;
                      cbf.d2fdxdy += p.v % d2edxdy;
                      cbf.d2fdxdz += p.v % d2edxdz;
                      cbf.d2fdydz += p.v % d2edydz;
                      cbf.d3fdxdydz += p.v % d3edxdydz;
                    }
                }
            }
        }
    }

  for ( size_t i = 0; i < spheres.size(); i++ )
    {
      sphere_t& s = spheres[i];
      vect3d r  = pos - s.pos;
      prec_t r2 = r.magnitude();
      prec_t rlen = prec_t_sqrt (r2);
      if ( r2 < s.rr )
        {
          epot.f += s.q_by_r;
          epot.f += 4.0 * PREC_PI * s.rho * (s.rr - r2) / 6.0;
          epot.dfdx += -4.0 * PREC_PI * s.rho * 2.0 * r.x / 6.0;
          epot.dfdy += -4.0 * PREC_PI * s.rho * 2.0 * r.y / 6.0;
          epot.dfdz += -4.0 * PREC_PI * s.rho * 2.0 * r.z / 6.0;
          prec_t f = s.q_by_r3;
          cef.f += f * r;
          cef.dfdx.x += f;
          cef.dfdy.y += f;
          cef.dfdz.z += f;
        }
      else
        {
          prec_t r3 = r2 * rlen;
          prec_t r5 = r3 * r2;
          prec_t r7 = r5 * r2;
          prec_t r9 = r7 * r2;

          epot.f += s.q / rlen;
          vect3d ef = (s.q / r3) * r;
          cef.f += ef;
          epot.dfdx += -ef.x;
          epot.dfdy += -ef.y;
          epot.dfdz += -ef.z;
          prec_t f = s.q / r5;
          vect3d dedx = f * (vect3d(r2, 0, 0) - 3 * r.x * r);
          vect3d dedy = f * (vect3d(0, r2, 0) - 3 * r.y * r);
          vect3d dedz = f * (vect3d(0, 0, r2) - 3 * r.z * r);
          cef.dfdx += dedx;
          cef.dfdy += dedy;
          cef.dfdz += dedz;
          epot.d2fdxdy += -dedy.x;
          epot.d2fdxdz += -dedy.x;
          epot.d2fdydz += -dedz.y;
          f = 3 * s.q / r7;
          vect3d d2edxdy = f * (5 * r.x * r.y * r - vect3d (r2 * r.y, r2 * r.x, 0));
          vect3d d2edxdz = f * (5 * r.x * r.z * r - vect3d (r2 * r.z, 0, r2 * r.x));
          vect3d d2edydz = f * (5 * r.y * r.z * r - vect3d (0, r2 * r.z, r2 * r.y));
          cef.d2fdxdy += d2edxdy;
          cef.d2fdxdz += d2edxdz;
          cef.d2fdydz += d2edydz;
          epot.d3fdxdydz += -d2edydz.x;
          f = 15 * s.q / r9;
          vect3d d3edxdydz = f * (vect3d(r.y * r.z * r2, r.x * r.z * r2, r.x * r.y * r2) - 7 * r.x * r.y * r.z * r);
          cef.d3fdxdydz += d3edxdydz;
        }
    }
  
  epot.mul (EMOD);
  cef.mul (EMOD);
  cbf.mul (BMOD);

#if 0
  cbf.dfdx.clear();
  cbf.dfdy.clear();
  cbf.dfdz.clear();
  cbf.d2fdxdy.clear();
  cbf.d2fdxdz.clear();
  cbf.d2fdydz.clear();
  cbf.d3fdxdydz.clear();
#endif

  out_epot.add (epot);
  out_cef.add (cef);
  out_cbf.add (cbf);
}

void PMap::addEField (const vect3d& pos, vect3d& ef)
{
  CubicScalar epot;
  CubicVector cbf, cef;
  cef.clear();
  addCubicValues (pos, epot, cbf, cef);
  ef += cef.f;
}

StaticPMap::StaticPMap (PMap& pmap, prec_t Emul) : pmap(pmap), Emul(Emul)
{
}

void StaticPMap::addFields (const vect3d& pos, vect3d& bfield, vect3d& efield)
{
  vect3d ef;
  ef.clear();
  pmap.addEField (pos, ef);
  efield += ef * Emul;
}

void StaticPMap::addPotential (const vect3d& pos, prec_t& potential)
{
  potential += Emul * pmap.getPotential (pos);
}

void StaticPMap::draw(Screen& screen, const Screen::color& c)
{
}

std::string StaticPMap::get_povray_script ()
{
  return std::string();
}

bool StaticPMap::inContact (const vect3d& pos)
{
  return pmap.inContact (pos);
}

StaticPMapCrudeEfield::StaticPMapCrudeEfield (PMap& pmap) : pmap(pmap)
{
}

void StaticPMapCrudeEfield::addFields (const vect3d& pos, vect3d& bfield, vect3d& efield)
{
  pmap.addCrudeEField (pos, efield);
}

void StaticPMapCrudeEfield::addPotential (const vect3d& pos, prec_t& potential)
{
  potential += pmap.getCrudePotential (pos);
}

void StaticPMapCrudeEfield::draw(Screen& screen, const Screen::color& c)
{
}

std::string StaticPMapCrudeEfield::get_povray_script ()
{
  return std::string();
}

bool StaticPMapCrudeEfield::inContact (const vect3d& pos)
{
  return pmap.inContact (pos);
}
